home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / FORTRAN1.LZH / KURV1.FOR < prev    next >
Text File  |  1988-02-08  |  6KB  |  181 lines
  1.       SUBROUTINE KURV1 ( N, X, Y, SLP1, SLPN, XP, YP, TEMP, S,
  2.      $ SIGMA )
  3. C*
  4. C*              *******************************************
  5. C*              *                                         *
  6. C*              *                   KURV1                 *
  7. C*              *                                         *
  8. C*              *******************************************
  9. C*
  10. C*      KURV1 - CLINE'S ROUTINE TO SET UP PIECEWISE CUBIC SPLINE UNDER
  11. C*           TENSION, FOR EVALUATION BY ROUTINE KURV2.
  12. C*
  13. C*      AUTHOR - A. K. CLINE
  14. C*               NATIONAL CENTER FOR ATMOSPHERIC RESEARCH
  15. C*               PO BOX 1470
  16. C*               BOULDER, COLORADO     80302
  17. C*
  18. C*      CODED BY - A. E. RAGOSTA        415-694-6235
  19. C*               TR18
  20. C*               AMES RSCH CTR
  21. C*               MOFFETT FIELD, CALIF  94035
  22. C*
  23. C*      INPUT ARGUMENTS
  24. C*           N     - NUMBER OF POINTS TO BE FIT
  25. C*           X     - ARRAY OF X VALUES
  26. C*           Y     - ARRAY OF Y VALUES
  27. C*           SLP1  - SLOPE AT FIRST POINT (DEGREES, COUNTER-CLOCKWISE FROM
  28. C*                       POSITIVE X AXIS)
  29. C*           SLPN  - SLOPE AT LAST POINT (DEGREES, ...)
  30. C*           SIGMA - TENSION FACTOR (IF THIS VALUE IS NEGATIVE, THE END
  31. C*                     POINT SLOPES WILL BE CALCULATED, IF POSITIVE, THEY
  32. C*                     SHOULD BE INPUT IN SLP1 AND SLP2.  A TYPICAL VALUE
  33. C*                     IS 1. )
  34. C*           TEMP  - SCRATCH WORK AREA
  35. C*
  36. C*      OUTPUT ARGUMENTS
  37. C*           XP    - CURVATURE PARAMETERS FOR KURV2
  38. C*           YP    - CURVATURE PARAMETERS FOR KURV2
  39. C*           S     - ARC LENGTH OF CURVE
  40. C*
  41. C*      COMMON BLOCKS REFERENCED :
  42. C*           NONE
  43. C*
  44. C*      FILES REFERENCED :
  45. C*           NONE
  46. C*
  47. C*     EXTERNAL SUBPROGRAMS REFERENCED :
  48. C*          SIN, SQRT, COS, ABS, FLOAT, EXP, ATAN2
  49. C*
  50. C*      VERSION I        29 DEC 1981
  51. C*
  52. C***********************************************************************
  53. C*
  54.       DIMENSION X(N), Y(N), YP(N), XP(N), TEMP(N)
  55. C
  56.       DEGRAD = 0.01745329
  57.       NM1   = N - 1
  58.       NP1   = N + 1
  59.       DELX1 = X(2) - X(1)
  60.       DELY1 = Y(2) - Y(1)
  61.       DELS1 = SQRT ( DELX1**2 + DELY1**2 )
  62.       DX1   = DELX1 / DELS1
  63.       DY1   = DELY1 / DELS1
  64. C
  65. C --- CALCULATE SLOPES IF REQUESTED
  66. C
  67.       IF ( SIGMA .LT. 0. ) GO TO 70
  68.       SLPP1 = SLP1 * DEGRAD
  69.       SLPPN = SLPN * DEGRAD
  70.    10 XP(1) = DX1 - COS ( SLPP1 )
  71. C
  72. C --- SET UP RIGHT HAND SIDE OF TRIDIAG
  73. C
  74.       YP(1) = DY1 - SIN ( SLPP1 )
  75.       TEMP(1) = DELS1
  76.       S     = DELS1
  77.       IF ( N .EQ. 2 ) GO TO 30
  78.       DO 20 I = 2, NM1
  79.          DELX2 = X(I+1) - X(I)
  80.          DELY2 = Y(I+1) - Y(I)
  81.          DELS2 = SQRT ( DELX2**2 + DELY2**2 )
  82.          DX2   = DELX2 / DELS2
  83.          DY2   = DELY2 / DELS2
  84.          XP(I) = DX2 - DX1
  85.          YP(I) = DY2 - DY1
  86.          TEMP(I) = DELS2
  87.          DELX1 = DELX2
  88.          DELY1 = DELY2
  89.          DELS1 = DELS2
  90.          DX1   = DX2
  91.          DY1   = DY2
  92.          S     = S + DELS1
  93.    20    CONTINUE
  94.    30 XP(N)  = COS ( SLPPN ) - DX1
  95.       YP(N)  = SIN ( SLPPN ) - DY1
  96. C
  97. C --- DENORMALIZE TENSION FACTOR
  98. C
  99.       SIGMAP = ABS ( SIGMA ) * FLOAT ( N-1 ) / S
  100. C
  101. C --- FORWARD ELIMINATION ON TRIDIAGONAL
  102. C
  103.       DELS   = SIGMAP * TEMP(1)
  104.       EXPS   = EXP ( DELS )
  105.       SINHS  = .5 * ( EXPS - 1./EXPS )
  106.       SINHIN = 1./( TEMP(1) * SINHS )
  107.       DIAG1  = SINHIN * ( DELS * .5 * ( EXPS + 1./EXPS ) - SINHS )
  108.       DIAGIN = 1./DIAG1
  109.       XP(1)  = DIAGIN * XP(1)
  110.       YP(1)  = DIAGIN * YP(1)
  111.       SPDIAG = SINHIN * ( SINHS - DELS )
  112.       TEMP(1) = DIAGIN * SPDIAG
  113.       IF ( N .EQ. 2 ) GO TO 50
  114.       DO 40 I = 2, NM1
  115.          DELS   = SIGMAP * TEMP(I)
  116.          EXPS   = EXP ( DELS )
  117.          SINHS  = .5 * ( EXPS - 1./EXPS )
  118.          SINHIN = 1./( TEMP(I) * SINHS )
  119.          DIAG2  = SINHIN * ( DELS * ( .5 * ( EXPS + 1./EXPS )) - SINHS )
  120.          DIAGIN = 1./( DIAG1 + DIAG2 - SPDIAG * TEMP(I-1))
  121.          XP(I)  = DIAGIN * ( XP(I) - SPDIAG * XP(I-1))
  122.          YP(I)  = DIAGIN * ( YP(I) - SPDIAG * YP(I-1))
  123.          SPDIAG = SINHIN * ( SINHS - DELS )
  124.          TEMP(I) = DIAGIN * SPDIAG
  125.          DIAG1  = DIAG2
  126.    40    CONTINUE
  127.    50 DIAGIN = 1./( DIAG1 - SPDIAG * TEMP(NM1))
  128.       XP(N)  = DIAGIN * ( XP(N) - SPDIAG * XP(NM1))
  129.       YP(N)  = DIAGIN * ( YP(N) - SPDIAG * YP(NM1))
  130. C
  131. C --- PERFORM SUBSTITUTIONS FOR COEFFICIENTS
  132. C
  133.       DO 60 I = 2, N
  134.          IBAK = NP1 - I
  135.          XP(IBAK) = XP(IBAK) - TEMP(IBAK) * XP(IBAK+1)
  136.          YP(IBAK) = YP(IBAK) - TEMP(IBAK) * YP(IBAK+1)
  137.    60    CONTINUE
  138.       GO TO 90
  139. C
  140. C --- SECOND ORDER INTERPOLATION FOR ENDPOINTS
  141. C      (IF NO SLOPES SPECIFIED)
  142. C
  143.    70 IF ( N .EQ. 2 ) GO TO 80
  144.       DELS2 = SQRT (( X(3) - X(2))**2 + ( Y(3) - Y(2))**2 )
  145.       DELS12 = DELS1 + DELS2
  146.       C1    = -(DELS12 + DELS1)/DELS12/DELS1
  147.       C2    = DELS12/DELS1/DELS2
  148.       C3    = -DELS1 / ( DELS12 * DELS2 )
  149.       SX    = C1 * X(1) + C2 * X(2)  +  C3 * X(3)
  150.       SY    = C1 * Y(1)  +  C2 * Y(2)  +  C3 * Y(3)
  151.       SLPP1 = ATAN2 ( SY, SX )
  152.       SLP1  = SLPP1 / DEGRAD +180.
  153.       DELNM1= SQRT (( X(N-2) - X(NM1))**2 + ( Y(N-2) - Y(NM1))**2 )
  154.       DELN  = SQRT (( X(NM1) - X(N))**2 + ( Y(NM1) - Y(N))**2 )
  155.       DELNN = DELNM1 + DELN
  156.       C1    = ( DELNN + DELN ) / ( DELNN * DELN )
  157.       C2    = -DELNN / ( DELN * DELNM1 )
  158.       C3    = DELN / ( DELNN * DELNM1 )
  159.       SX    = C3 * X(N-2)  +  C2 * X(NM1)  +  C1 * X(N)
  160.       SY    = C3 * Y(N-2)  +  C2 * Y(NM1)  +  C1 * Y(N)
  161.       SLPPN = ATAN2 ( SY, SX )
  162.       SLPN  = SLPPN / DEGRAD
  163.       IF ( SLPN .LT. 0. )SLPN = SLPN + 360.
  164.       GO TO 10
  165. C
  166. C --- TWO POINTS ONLY, RETURN A STRAIGHT LINE
  167. C
  168.    80 XP(1) = 0.
  169.       XP(2) = 0.
  170.       YP(1) = 0.
  171.       YP(2) = 0.
  172.       SLP1 = ATAN2 ((Y(2)-Y(1)),(X(2)-X(1))) / DEGRAD
  173.       SLPN = SLP1
  174.       IF ( SLPN .LT. 0. ) SLPN = SLPN + 360.
  175.       SLP1 = SLP1 + 180.
  176.    90 RETURN
  177.       END
  178. C
  179. C---END KURV1
  180. C
  181.